Compressive Wave Computation

This paper considers large-scale simulations of wave propagation phenomena. We argue that it is possible to accurately compute a wavefield by decomposing it onto a largely incomplete set of eigenfunctions of the Helmholtz operator, chosen at random, and that this provides a natural way of parallelizing wave simulations for memory-intensive applications. Where a standard eigenfunction expansion in general fails to be accurate if a single term is missing, a sparsity-promoting `1 minimization problem can vastly enhance the quality of synthesis of a wavefield from low-dimensional spectral information. This phenomenon may be seen as “compressive sampling in the Helmholtz domain”, and has recently been observed to have a bearing on the performance of data extrapolation techniques in seismic imaging [41]. This paper shows that `1-Helmholtz recovery also makes sense for wave computation, and identifies a regime in which it is provably effective: the one-dimensional wave equation with coefficients of small bounded variation. Under suitable assumptions we show that the number of eigenfunctions needed to evolve a sparse wavefield defined on N points, accurately with very high probability, is bounded by C(η) · logN · log logN, where C(η) is related to the desired accuracy η and can be made to grow at a much slower rate than N when the solution is sparse. The PDE estimates that underlie this result are new to the authors’ knowledge and may be of independent mathematical interest; they include an L estimate for the wave equation, an L∞ − L estimate of extension of eigenfunctions, and a bound for eigenvalue gaps in Sturm-Liouville problems. In practice, the compressive strategy makes sense because the computation of eigenfunctions can be assigned to different nodes of a cluster in an embarrassingly parallel way. Numerical examples are presented in one spatial dimension and show that as few as 10 percents of all eigenfunctions can suffice for accurate results. Availability of a good preconditioner for the Helmholtz equation is important and also discussed in the paper. Finally, we argue that the compressive viewpoint suggests a competitive parallel algorithm for an adjoint-state inversion method in reflection seismology. Acknowledgements. We would like to thank Ralph Smith and Jim Berger of the Statistical and Applied Mathematical Sciences Institute (SAMSI) for giving us the opportunity to visit the Institute, which catalyzed the completion of this project. We are indebted to Felix Herrmann for early discussions, Paul Rubin for help on a probability question related to Proposition 4, and Lexing Ying for help with a point in the numerical implementation of our preconditioner. Thanks are also due to Emmanuel Candès, Mark Embree, Jalal Fadili, Josselin Garnier, Mauro Maggioni, Justin Romberg, and William Symes for useful discussions. L.D. is supported in part by a grant from the National Science Foundation.